package soeasy.util;

import java.awt.Color;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.util.List;
import java.util.StringTokenizer;
import java.util.Vector;

import org.jfree.data.Range;
import org.jfree.data.xy.XYDataItem;
import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;

import soeasy.model.SOEASYParameters;
import umontreal.iro.lecuyer.charts.XYLineChart;
import umontreal.iro.lecuyer.charts.XYListSeriesCollection;
import cern.colt.list.DoubleArrayList;

public class CuSum {

	private Range offsetRange = null;

	// general properties
	private double prestimulusMean;
	private double errorBoxUpperBound;
	private double errorBoxLowerBound;
	private List<Event> eventList;
	private String diagramType;
	private String units;

	/**
	 * The size of each bin.
	 */
	protected double binSize = 0.5;

	public DigitalizedXYSeries cusumXYSeries;

	public DigitalizedXYSeries derivativeXYSeries;

	private double derivativePrestimulusMean = 0.0;

	private double derivativeStandardDeviation = 0.0;

	private double derivativeUpperControlLimit = 0.0;

	private double derivativeLowerControlLimit = 0.0;

	public DigitalizedXYSeries pureDerivativeXYSeries;

	private double pureDerivativePrestimulusMean = 0.0;

	private double pureDerivativeStandardDeviation = 0.0;

	private double pureDerivativeUpperControlLimit = 0.0;

	private double pureDerivativeLowerControlLimit = 0.0;

	private int smooth = 1;

	private double movingAverageInterval = 4.0;

	/**
	 * Creates an empty CuSum object
	 */
	public CuSum() {
		cusumXYSeries = new DigitalizedXYSeries(this.offsetRange, 0.5);
		diagramType = "x";
		units = "x";
	}

	public CuSum(PeriStimulusTimeHistogram psth, int smooth,
			double movingAverageInterval) {
		this.offsetRange = psth.getOffsetRange();
		this.smooth = smooth;
		this.movingAverageInterval = movingAverageInterval;
		//
		// this.binSize = 1.0;
		cusumXYSeries = new DigitalizedXYSeries(this.offsetRange, this.binSize);
		derivativeXYSeries = new DigitalizedXYSeries(this.offsetRange,
				this.binSize);
		pureDerivativeXYSeries = new DigitalizedXYSeries(this.offsetRange,
				this.binSize);
		diagramType = "PSTH";
		units = "counts.ms";

		// calculate the pre-stimulus mean
		this.prestimulusMean = psth.getPrestimulusMean();

		// calculate the error box (lower and upper bounds)
		errorBoxUpperBound = this.prestimulusMean;
		errorBoxLowerBound = this.prestimulusMean;

		// the first cusum value is 0.0
		double currentCUSUM = 0.0;
		cusumXYSeries.add(offsetRange.getLowerBound(), currentCUSUM);
		currentCUSUM += calculateDifference(psth, 0);

		int endIndex = psth.getBins().size() - 1;
		for (int index = 1; index <= endIndex; index++) {
			currentCUSUM += calculateDifference(psth, index);

			double xValue = offsetRange.getLowerBound()
					+ (index * psth.getBinSize());
			cusumXYSeries.add(xValue, currentCUSUM);

			// find the highest bin
			if ((currentCUSUM > errorBoxUpperBound)
					&& (index <= psth.getZeroIndex())) {
				errorBoxUpperBound = currentCUSUM;
			}

			// and find the lowest bin
			if ((currentCUSUM < errorBoxLowerBound)
					&& (index <= psth.getZeroIndex())) {
				errorBoxLowerBound = currentCUSUM;
			}
		}

		double upperDifference = Math.abs(errorBoxUpperBound - prestimulusMean);
		double lowerDifference = Math.abs(errorBoxLowerBound - prestimulusMean);
		if (upperDifference > lowerDifference) {
			errorBoxLowerBound = prestimulusMean - upperDifference;
		} else {
			errorBoxUpperBound = prestimulusMean + lowerDifference;
		}

		calculateDerivative();

		calculatePureDerivativePrestimulusMean();
		calculateDerivativePrestimulusMean();

		calculatePureDerivativeStandardDeviation();
		calculateDerivativeStandardDeviation();

		calculatePureDerivativeControlLimits();
		calculateDerivativeControlLimits();
	}

	public CuSum(PeristimulusFrequencygram psf, int smooth,
			double movingAverageInterval) {
		this.offsetRange = psf.offsetRange;
		this.binSize = psf.getBinSize();
		this.smooth = smooth;
		this.movingAverageInterval = movingAverageInterval;

		// initialization
		cusumXYSeries = new DigitalizedXYSeries(this.offsetRange, this.binSize);
		derivativeXYSeries = new DigitalizedXYSeries(this.offsetRange,
				this.binSize);
		pureDerivativeXYSeries = new DigitalizedXYSeries(this.offsetRange,
				this.binSize);
		diagramType = "PSF";
		units = "Hz.ms";

		calculatePSFCUSUM(psf);

		calculateDerivative();

		calculatePureDerivativePrestimulusMean();
		calculateDerivativePrestimulusMean();

		calculatePureDerivativeStandardDeviation();
		calculateDerivativeStandardDeviation();

		calculatePureDerivativeControlLimits();
		calculateDerivativeControlLimits();
	}

	private void calculatePSFCUSUM(PeristimulusFrequencygram psf) {
		DigitalizedXYSeries dischargeRates = psf.getDischargeRates();

		// set the pre-stimulus mean
		this.prestimulusMean = psf.getMeanPrestimulusDischargeRate();
		//
		// calculate the error box (lower and upper bounds)
		errorBoxUpperBound = this.prestimulusMean;
		errorBoxLowerBound = this.prestimulusMean;

		// the first cusum value is set to 0.0
		double psfCUSUM = 0.0;
		cusumXYSeries.add(offsetRange.getLowerBound(), psfCUSUM);

		// then from the second bin to the end of the bins calculate cusums
		for (double xValue = offsetRange.getLowerBound() + binSize; xValue <= offsetRange
				.getUpperBound(); xValue += binSize) {
			psfCUSUM += calculateDifference(dischargeRates, xValue);
			cusumXYSeries.add(xValue, psfCUSUM);
			//
			// find the highest bin
			if ((psfCUSUM > errorBoxUpperBound) && (xValue <= 0.0)) {
				errorBoxUpperBound = psfCUSUM;
			}

			// and find the lowest bin
			if ((psfCUSUM < errorBoxLowerBound) && (xValue <= 0.0)) {
				errorBoxLowerBound = psfCUSUM;
			}
		}
	}

	// derivative(y) = dy/dx
	private void calculateDerivative() {
		if (cusumXYSeries.getItemCount() > 1) {
			// calculate derivatives
			for (double xValue = offsetRange.getLowerBound() + binSize; xValue <= offsetRange
					.getUpperBound(); xValue += binSize) {
				try {
					double x1 = xValue - binSize;
					double x2 = xValue;
					double y1 = cusumXYSeries.getItemsAt(x1).get(0).getYValue();
					double y2 = cusumXYSeries.getItemsAt(x2).get(0).getYValue();

					double derivative = (y2 - y1) / 2; // (x2 - x1);
					derivativeXYSeries.add(x2, derivative);
					pureDerivativeXYSeries.add(x2, derivative);
				} catch (Exception e) {
					// System.out.println("xValue: " + xValue);
				}
			}

			derivativeXYSeries = calculateMovingAvarage(derivativeXYSeries,
					smooth);
		}
	}

	private void calculatePureDerivativePrestimulusMean() {
		Averager averager = new Averager();
		for (double xValue = offsetRange.getLowerBound() + binSize; xValue <= 0.0; xValue += binSize) {
			List<XYDataItem> itemsAt = pureDerivativeXYSeries
					.getItemsAt(xValue);
			double yValue = itemsAt.get(0).getYValue();
			averager.update(yValue);
		}
		this.pureDerivativePrestimulusMean = averager.getAverage();
	}

	private void calculateDerivativePrestimulusMean() {
		Averager averager = new Averager();
		for (double xValue = offsetRange.getLowerBound() + binSize; xValue <= 0.0; xValue += binSize) {
			List<XYDataItem> itemsAt = derivativeXYSeries.getItemsAt(xValue);
			double yValue = itemsAt.get(0).getYValue();
			averager.update(yValue);
		}
		this.derivativePrestimulusMean = averager.getAverage();
	}

	private void calculatePureDerivativeStandardDeviation() {
		Averager squareAverager = new Averager();
		for (double xValue = offsetRange.getLowerBound() + binSize; xValue <= 0.0; xValue += binSize) {
			List<XYDataItem> itemsAt = pureDerivativeXYSeries
					.getItemsAt(xValue);
			double yValue = itemsAt.get(0).getYValue();
			double diff = yValue - this.pureDerivativePrestimulusMean;
			squareAverager.update(diff * diff);
		}
		this.pureDerivativeStandardDeviation = Math.sqrt(squareAverager
				.getAverage());
	}

	private void calculateDerivativeStandardDeviation() {
		Averager squareAverager = new Averager();
		for (double xValue = offsetRange.getLowerBound() + binSize; xValue <= 0.0; xValue += binSize) {
			List<XYDataItem> itemsAt = derivativeXYSeries.getItemsAt(xValue);
			double yValue = itemsAt.get(0).getYValue();
			double diff = yValue - this.derivativePrestimulusMean;
			squareAverager.update(diff * diff);
		}
		this.derivativeStandardDeviation = Math.sqrt(squareAverager
				.getAverage());
	}

	private void calculatePureDerivativeControlLimits() {
		this.pureDerivativeUpperControlLimit = pureDerivativePrestimulusMean
				+ pureDerivativeStandardDeviation * 3;
		this.pureDerivativeLowerControlLimit = pureDerivativePrestimulusMean
				- pureDerivativeStandardDeviation * 3;
	}

	private void calculateDerivativeControlLimits() {
		this.derivativeUpperControlLimit = derivativePrestimulusMean
				+ derivativeStandardDeviation * 3;
		this.derivativeLowerControlLimit = derivativePrestimulusMean
				- derivativeStandardDeviation * 3;
	}

	public Range getDerivationInterval(double latency) {
		double derivation = derivativeXYSeries.getItemsAt(latency).get(0)
				.getYValue();

		double interval = SOEASYParameters.getInstance().getParameter(
				SOEASYParameters.DERIVATIVE_FAULT_TOLERANCE);
		double min = derivation - (interval / 2);
		double max = derivation + (interval / 2);

		return new Range(min, max);
	}

	public double getDerivativePrestimulusMean() {
		return derivativePrestimulusMean;
	}

	private DigitalizedXYSeries calculateMovingAvarage(
			DigitalizedXYSeries inputs, int smooth) {
		DigitalizedXYSeries averageValues = new DigitalizedXYSeries(
				offsetRange, binSize);

		for (double time = offsetRange.getLowerBound(); time <= offsetRange
				.getUpperBound(); time += binSize) {
			double fromXValue = time - movingAverageInterval / 2;
			double toXValue = time + movingAverageInterval / 2;

			double averageDischargeRate = inputs.getAverageYValueBetween(
					fromXValue, toXValue);
			if (!Double.isNaN(averageDischargeRate)) {
				averageValues.add(time, averageDischargeRate);
			}
		}

		// smooth
		if (smooth > 0) {
			averageValues = calculateMovingAvarage(averageValues, (smooth - 1));
		}

		//
		return averageValues;
	}

	/**
	 * Calculates and returns the total difference of all the values at xValue
	 * from the pre-stimulus mean.
	 */
	private double calculateDifference(DigitalizedXYSeries dischargeRates,
			double xValue) {
		// get the averager for "bin xValue"
		Averager averagerAtX = dischargeRates.getAveragerAt(xValue);

		// get number of values at "bin xValue"
		double elementCount = averagerAtX.getElementCount();

		// get the average discharge rate value at "bin xValue"
		double averageDischargeRateValue = averagerAtX.getAverage();
		double totalDiff = (averageDischargeRateValue - prestimulusMean)
				* elementCount;
		return totalDiff;
	}

	private double calculateDifference(PeriStimulusTimeHistogram psth, int index) {
		int currentBinValue = psth.getBinValue(index);
		double totalDiff = currentBinValue - prestimulusMean;
		return totalDiff;
	}

	public double getCuSumAt(double time) {
		double cusum = Double.NaN;
		List<XYDataItem> itemsAtTime = cusumXYSeries.getItemsAt(time);
		if (!itemsAtTime.isEmpty()) {
			cusum = itemsAtTime.get(0).getYValue();
		}
		return cusum;
	}

	public List<Event> getEvents() {
		if (this.eventList == null) {
			this.eventList = new Vector<Event>();
		}
		return this.eventList;
	}

	public double getPrestimulusMean() {
		return prestimulusMean;
	}

	public double getErrorBoxUpperBound() {
		return errorBoxUpperBound;
	}

	public double getErrorBoxLowerBound() {
		return errorBoxLowerBound;
	}

	public double getDerivativeStandardDeviation() {
		return derivativeStandardDeviation;
	}

	public double getDerivativeUpperControlLimit() {
		return derivativeUpperControlLimit;
	}

	public double getDerivativeLowerControlLimit() {
		return derivativeLowerControlLimit;
	}

	public void toFile(String fileName) {
		try {
			FileWriter fileWriter = new FileWriter(fileName);
			BufferedWriter bufferedWriter = new BufferedWriter(fileWriter);
			bufferedWriter.write("" + prestimulusMean + "\n");
			bufferedWriter.write("" + errorBoxLowerBound + "\n");
			bufferedWriter.write("" + errorBoxUpperBound + "\n");
			for (int index = 0; index < cusumXYSeries.getItemCount(); index++) {
				XYDataItem dataItem = cusumXYSeries.getDataItem(index);
				bufferedWriter.write(dataItem.getXValue() + " "
						+ dataItem.getYValue() + "\n");

			}
			bufferedWriter.close();
		} catch (Exception e) {
			e.printStackTrace();
		}
	}

	public void readFile(String fileName) {
		try {
			FileReader fileReader = new FileReader(fileName);
			BufferedReader bufferedReader = new BufferedReader(fileReader);
			//
			String line = bufferedReader.readLine();
			prestimulusMean = Double.parseDouble(line);
			//
			line = bufferedReader.readLine();
			errorBoxLowerBound = Double.parseDouble(line);
			//
			line = bufferedReader.readLine();
			errorBoxUpperBound = Double.parseDouble(line);
			//
			while ((line = bufferedReader.readLine()) != null) {
				StringTokenizer tokenizer = new StringTokenizer(line);
				double xValue = Double.parseDouble(tokenizer.nextToken());
				double yValue = Double.parseDouble(tokenizer.nextToken());
				cusumXYSeries.add(xValue, yValue);
			}
			//
			bufferedReader.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	public XYLineChart getXYChart() {
		return getXYChart(offsetRange);
	}

	public XYLineChart getXYChart(Range offsetRange) {
		Range cusumRange = getCUSUMRange(offsetRange);
		return getXYChart(offsetRange, cusumRange);
	}

	public XYLineChart getXYChart(Range offsetRange, Range cusumRange) {
		XYSeriesCollection xySeriesCollection = new XYSeriesCollection();

		xySeriesCollection.addSeries(this.cusumXYSeries.toXYSeries(""));
		xySeriesCollection.addSeries(getYSeries(errorBoxUpperBound,
				"Error Box upper bound"));
		xySeriesCollection.addSeries(getYSeries(errorBoxLowerBound,
				"Error Box lower bound"));
		xySeriesCollection.addSeries(getYSeries(prestimulusMean,
				"Prestimulus mean"));

		XYLineChart chart = new XYLineChart(diagramType + "-CUSUM",
				"Time (ms)", diagramType + "-CUSUM (" + units + ")",
				xySeriesCollection);

		// Customizes the data plot
		XYListSeriesCollection collec = chart.getSeriesCollection();
		collec.setColor(0, Color.RED); // cusum
		collec.setColor(1, Color.BLUE); // error box upper bound
		collec.setColor(2, Color.BLUE); // error box lower bound
		collec.setColor(3, Color.GRAY); // prestimulus mean
		collec.setColor(4, Color.WHITE); // stimulus
		// now set white for every event

		if (cusumRange != null) {
			double[] bounds = { offsetRange.getLowerBound(),
					offsetRange.getUpperBound(),
					cusumRange.getLowerBound() - 150,
					cusumRange.getUpperBound() + 150 };
			chart.setManualRange(bounds);
		} else {
			chart.setAutoRange();
		}

		return chart;
	}

	protected Range getCUSUMRange(Range offsetRange) {
		double minimumCUSUM = Double.MAX_VALUE;
		double maximumCUSUM = Double.MIN_VALUE;

		double from = offsetRange.getLowerBound();
		double to = offsetRange.getUpperBound();
		for (double xValue = from; xValue <= to; xValue += 0.5) {
			List<XYDataItem> items = cusumXYSeries.getItemsAt(xValue);
			if (items.size() > 0) {
				XYDataItem item = items.get(0);
				double cusum = item.getYValue();
				if (cusum < minimumCUSUM) {
					minimumCUSUM = cusum;
				} else if (cusum > maximumCUSUM) {
					maximumCUSUM = cusum;
				}

			}
		}
		return new Range(minimumCUSUM, maximumCUSUM);
	}

	public XYLineChart getXYChartAutoRange() {
		XYSeriesCollection xySeriesCollection = new XYSeriesCollection();

		xySeriesCollection.addSeries(this.cusumXYSeries.toXYSeries(""));
		xySeriesCollection.addSeries(getYSeries(errorBoxUpperBound,
				"Error Box upper bound"));
		xySeriesCollection.addSeries(getYSeries(errorBoxLowerBound,
				"Error Box lower bound"));
		xySeriesCollection.addSeries(getYSeries(prestimulusMean,
				"Prestimulus mean"));

		// stimulus
		xySeriesCollection.addSeries(getXSeries(0));

		XYLineChart chart = new XYLineChart(diagramType + "-CUSUM",
				"Time (ms)", diagramType + "-CUSUM (" + units + ")",
				xySeriesCollection);

		// Customizes the data plot
		XYListSeriesCollection collec = chart.getSeriesCollection();
		collec.setColor(0, Color.RED); // cusum
		collec.setColor(1, Color.BLUE); // error box upper bound
		collec.setColor(2, Color.BLUE); // error box lower bound
		collec.setColor(3, Color.GRAY); // prestimulus mean
		collec.setColor(4, Color.WHITE); // stimulus

		chart.setAutoRange();

		return chart;
	}

	private XYSeries getXSeries(double xValue) {
		XYSeries xySeries = new XYSeries("");
		for (int y = -400; y < 400; y++) {
			xySeries.add(xValue, y);
		}
		return xySeries;
	}

	private XYSeries getYSeries(double yValue, String key) {
		XYSeries xySeries = new XYSeries(key);
		for (double x = offsetRange.getLowerBound(); x < offsetRange
				.getUpperBound(); x++) {
			xySeries.add(x, yValue);
		}
		return xySeries;
	}

	public double[] getYValues(double from, double to, double sensitivity) {
		DoubleArrayList yValuesList = new DoubleArrayList();
		for (double index = from; index <= to; index = index + sensitivity) {
			yValuesList.add(getCuSumAt(index));
		}
		return yValuesList.elements();
	}

	@Override
	public String toString() {
		StringBuffer buffer = new StringBuffer();
		buffer.append("[");
		buffer.append("error box lower: " + this.errorBoxLowerBound + "; ");
		buffer.append("prestimulus mean: " + this.prestimulusMean + "; ");
		buffer.append("error box upper: " + this.errorBoxUpperBound);
		buffer.append("]");
		return buffer.toString();
	}

	public DigitalizedXYSeries getCusumXYSeries() {
		return this.cusumXYSeries;
	}

	public DigitalizedXYSeries getDerivativeXYSeries() {
		return this.derivativeXYSeries;
	}

	public double getPureDerivativeAt(double latency) {
		return this.pureDerivativeXYSeries.getItemsAt(latency).get(0)
				.getYValue();
	}

	public double getDerivativeAt(double latency) {
		return this.derivativeXYSeries.getItemsAt(latency).get(0).getYValue();
	}

	public double getDerivativeAsStdDevAt(double latency) {
		double derivative = this.getDerivativeAt(latency);
		return (derivative - derivativePrestimulusMean)
				/ derivativeStandardDeviation;
	}

	/**
	 * mean squared error
	 * 
	 * @param simulatedPSFCUSUM
	 * @return
	 */
	public double calculatePrestimulsDerivationError(CuSum otherCUSUM) {
		Averager averager = new Averager();
		for (double xValue = offsetRange.getLowerBound(); xValue < 0.0; xValue += binSize) {
			double diff = this.getDerivativeAt(xValue)
					- otherCUSUM.getDerivativeAt(xValue);
			averager.update(diff * diff);
		}
		return averager.getAverage();
	}

	public boolean isPureBetweenControlLimits(double latency) {
		return !isPureBelowControlLimits(latency)
				&& !isPureAboveControlLimits(latency);
	}

	public boolean isBetweenControlLimits(double latency) {
		return !isBelowControlLimits(latency) && !isAboveControlLimits(latency);
	}

	public boolean isPureAboveControlLimits(double latency) {
		double derivative = getPureDerivativeAt(latency);
		return derivative > pureDerivativeUpperControlLimit;
	}

	private boolean isAboveControlLimits(double latency) {
		double derivative = getDerivativeAt(latency);
		return derivative > derivativeUpperControlLimit;
	}

	public boolean isPureBelowControlLimits(double latency) {
		double derivative = getPureDerivativeAt(latency);
		return derivative < pureDerivativeLowerControlLimit;
	}

	private boolean isBelowControlLimits(double latency) {
		double derivative = getDerivativeAt(latency);
		return derivative < derivativeLowerControlLimit;
	}

	public int compareAt(CuSum otherCUSUM, double latency) {
		int comparison = 0;

		boolean betweenControlLimits = this.isBetweenControlLimits(latency);
		if (betweenControlLimits) {
			if (otherCUSUM.isBetweenControlLimits(latency)) {
				comparison = 0;
			} else if (otherCUSUM.isAboveControlLimits(latency)) {
				comparison = 1;
			} else if (otherCUSUM.isBelowControlLimits(latency)) {
				comparison = -1;
			}
		} else {
			double diff = otherCUSUM.getDerivativeAsStdDevAt(latency)
					- this.getDerivativeAsStdDevAt(latency);
			double tolerance = 0.5;

			if (diff > tolerance) {
				comparison = 1;
			} else if (diff < -tolerance) {
				comparison = -1;
			} else {
				comparison = 0;
			}
		}

		return comparison;
	}
}
